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Abstract 

The Cellular Potts Model {CPM) is a robust, cell-level methodology for simulation of biological 
tissues and morphogenesis. Both tissue physiology and morphogenesis depend on diffusion of chem- 
ical morphogens in the extra-cellular fluid or matrix {ECM). Standard diffusion solvers applied to 
the cellular potts model use finite difference methods on the underlying CPM lattice. However, 
these methods produce a diffusing field tied to the underlying lattice, which is inaccurate in many 
biological situations in which cell or ECM movement causes advection rapid compared to diffusion. 
Finite difference schemes suffer numerical instabilities solving the resulting advection-diffusion 
equations. To circumvent these problems we simulate advection-diffusion within the framework of 
the CPM using off-lattice finite-difference methods. We define a set of generalized fluid particles 
which detach advection and diffusion from the lattice. Diffusion occurs between neighboring fluid 
particles by local averaging rules which approximate the Laplacian. Directed spin flips in the CPM 
handle the advective movement of the fluid particles. A constraint on relative velocities in the 
fluid explicitly accounts for fluid viscosity. We use the CPM to solve various diffusion examples 
including multiple instantaneous sources, continuous sources, moving sources and different bound- 
ary geometries and conditions to validate our approximation against analytical and established 
numerical solutions. We also verify the CPM results for Poiseuille flow and Taylor-Aris dispersion. 
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I. INTRODUCTION 



Advection-Diffusion equations (ADE) describe a broad range of natural phenomena. 

They occur in diverse field including physics Q , chemistry [2] , biology, geology ^, |^ and 

even in migration and epidemiology 5]. They describe the flow (deterministic) and the 

spread (stochastic) of a density (of a chemical, heat, charge) which a fluid or deformable 

solid carries. The simplest ADE is : 

c)ti — ^ 

— = DV 2 n - ~v ■ Vn, (1) 

dt w 

where n is the density of the transported substance, D its diffusion constant (here assumed 
uniform in space) and is the velocity field. The velocity field in turn couples to the 
pressure field of the medium through the Navier-Stokes equations. Though we can solve the 
problem analytically in steady state with simple boundary conditions, most physically rel- 
evant ADEs appear within sets of nonlinear coupled equations or with nontrivial boundary 
conditions where analytical solutions are not possible 6]. Hence a vast literature exists on 
low to solve ADEs. Most solvers use either finite difference (FD) or finite element (FE) 
3] schemes. Besides these deterministic approaches, several schemes use Lattice-Boltzmann 
(LB) methods like Flekkoy's method .9], Dawson's method 3], the moment-propagation 



method 



ADEs in general are difficult to solve in the absence of separation of diffu- 
sion and advection time scales or in the presence of moving boundaries. Most lattice-based 
methods locally refine the grid during solution to avoid instabilities. Moreover, explicit LB 
methods require time averaging o tt heto rq ue to avoid —ties Q. H^ce the eo.puta- 
tional cost of both LB and deterministic methods shoots up. Non-staggered FD grids may 
show grid-decoupling instabilities Q|. Also, all explicit methods require consideration of the 
general stability constraints from linear analysis, most notably the Neumann diffusive crite- 
rion linking the time step and the square of the grid size. In this article we try to address the 
problems associated with incorporating advection-diffusion in biologically-motivated, mul- 
tiscale simulations, specifically those which use the Cellular Potts Model (CPM) to model 
cell behaviors. 

Diffusion of morphogens and flow of extra-cellular matrix (ECM) are crucial to many 
biological phenomena, including wound healing, morphogenesis, e.g., during mesenchymal 



condensation or gastrulation 



131 ] and the immune response where cells emerge from the mi- 



crovasculature and migrate towards sites of inflammation to kill bacteria, other pathogens 



and cancer. The generic mechanisms common to all these processes are changes in cell veloc- 
ities (chemotaxis) or/and differentiation in response to the temporal and spatial variations 
of chemical morphogens. Other classic examples of diffusion-driven morphogenesis involve 
Turing instabilities. Turing instabilities arise due to different diffusion rates of two or more 
reacting chemicals resulting in competition between activation by a slow- diffusing chemi- 
cal (activator) and inhibition by a faster chemical (inhibitor). Pattern formation during 
morphogenesis due to Turing instabilities is a subject by itself 

Besides chemotaxis, the formation and rate of extension of pseudopods which crucially 
depends convective mass transport 14 , 3] also influences cell motility. Hydrodynamic shear 
can also increase cell-cell adhesion efficiency by increasing the number of binding receptors. 
Shear has a profound effect on neutrophil-platelet adhesion and neutrophil aggregation, key 
events in acute coronary syndromes like arterial thrombosis (fl^). Extensive work has shown 
that both fluid shear amplitude and shear exposure time modulate the interactions between 
polymorphonuclear leukocytes and colon carcinoma cells Gene expression and protein 
synthesis in endothelial cells also change upon application of arterial shear stresses In 
a prominent example, fluid shear allows optimal L-selectin-mediated leukocyte rolling only 
above a minimum threshold shear rate |l9j |. Hence multicellular modeling tools have to 
properly account for the advection-diffusion: diffusion influencing the spatial and temporal 
distribution of chemical morphogens, advection controlling the rates of cell collisions, de- 
formation, receptor-ligand bond formation [2]], adhesion and enhanced mixing of chemical 
morphogens. 

Simulations of the development of multicellular organisms take diverse mathematical ap- 



proaches: continuum FE-based models 



22j of react ion- diffusion which consider cell density 



as a continuous variable |23j, hybrid models like E-cell |2J] and cellular automaton ap- 
preaches Q. G.a Z ie r and Grader deve.oped the eelMevel CPM, an extension of the energy- 
based large-Q Potts model, for organogenesis simulations [26]. The basic CPM explains how 
surface binding energies drive cell movement and models cell sorting from an initial random 
distribution into different patterns depen ding on the cell adhesion coefficients at homotypic, 
heterotypic and cell-medium boundaries 271 ] . It also provides a platform on which to build 
simulations of a wide range of biological experiments by including additional mechanisms 
like directed active movement due to external fields, e.g. chemotaxis to a chemical field 
gradient, gravity or cell polarity. CPM applications include modeling mesenchymal conden- 
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sation 
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i2ai3o(, te life-cycle of Dictyostelium discoideum |35[, tumor growth 

vascular development [3^, immune response and limb growth j^J]. Unlike the simple 
Turing mechanism where cells have no feedback on the chemical field, most CPM implemen- 
tations include this feedback which can give rise to completely different patterning from the 
Turing mechanism. In the CPM, patterns can arise under the influence of a single chemi- 
cal field [13[ due to cell movement, biased by gradients in cell-cell adhesion and cell-ECM 
binding, which is impossible in Turing mechanism. This unique mechanism also differs from 
chemotaxis, which requires long-range cell movement j^ . 

All existing CPM implementations suffer from four main limitations: (1) They do 
not include viscous dissipation explicitly. Instead dissipation arises from the Metropolis- 
Boltzmann energy-minimization dynamics. This implicit dissipation makes viscosity hard 
to calculate or control. (2) They do not explicitly describe force transduction through cells, 
which arises through the volume constraint and surface constraint (if used) of the Hamil- 
tonian. (3) Aggregates of cells modeled with an ordinary CPM Hamiltonian are highly 
overdamped. Modeling the ECM or fluid as an array of generalized cells using the normal 
CPM produces a flow resembling the overdamped flow of biological fluid but this fluid slips 
at solid surfaces and exerts no shear force unlike a normal fluid. An alternative approach 
which describes the fluid as a single, large, unconstrained generalized cell produces nonlocal 
movement. Moreover, it cannot advect chemicals or transmit shear forces. (4) The CPM 
has no intrinsic concept of rigid-body motion. We will describe an algorithmic solution to 
this last problem, which makes the CPM look more like a FE simulation, in a future paper. 

All of these problems result from the CPM spins being tied to the underlying lattice, 
rather than to objects they describe. The solution in each case is to adopt an FE approach 
suited to the CPM, which takes the behavior off-lattice. Applying standard off-lattice meth- 
ods in the CPM has inherent problems. Since in the CPM cell movement occurs by boundary 
fluctuations, connecting normal FE fluid-solvers to the CPM at surfaces requires local grid 
refinement during cell movement to avoid numerical instabilities. Hence its computational 
cost is high. 

To introduce advection-diffusion in the CPM correctly and efficiently, we propose an off- 
lattice scheme consistent and in harmony with the CPM algorithm. We subdivide the ECM 
into small fluid particles having the normal properties of generalized cells like differential 
adhesivity, a surface constraint and a volume constraint. The fluid particles carry chemical 
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morphogens and we assume the concentrations are uniform inside them. Diffusion occurs 
between neighboring fluid particles by local averaging rules which approximate the Lapla- 
cian. Spin flips in the CPM occur only at the boundaries of the cells or particles. A force on 
the fluid in any direction due to pressure gradients or external forces biases the probability 
of spin flips and generates directed motion 3(|. We introduce viscosity into the CPM by 
explicitly in the Hamiltonian including a relative-velocity constraint between the neighbor- 
ing fluid particles. This scheme allows us to solve the ADEs and generates the creeping flow 
of a highly viscous fluid. 

In most biologically relevant regimes (e.g., E. Coli in water), we encounter low Reynolds 
number (Re) flow ~ 10~ 5 , with the typical diffusion coefficient of chemical morphogens being 
~ \0~ A [im 2 /s. Thus the Peclet number (defined as where R,v,D are typical size of the 
system, the typical fluid velocity and diffusion coefficient respectively) is as low as ~ f0~ 2 
(37 1 . We explicitly exclude blood flow in the circulatory system , where Re numbers (hence 
inertia) can be quite large and where specialized solution techniques already exist [? ] . Since 
most CPM simulations do not demand high precision, sophisticated methods like mimetic 
finite difference j^J become a computational bottle-neck without many advantages. On 
the other hand, our ADE scheme is very stable, with spatial resolution equal to the mean 
diameter of the fluid particles, which are much smaller than the simulated biological cells. 
Our scheme seamlessly integrates into the main Monte-Carlo loop of the CPM simulation 
and boundary conditions like absorbing boundaries or no-flux boundaries are simple to 
implement. More-over our off-lattice scheme does not require re-meshing. We discuss these 
issues in detail below. 

This paper mainly focuses on the validity and utility of the CPM ADE solver. Section 
2 briefly describes the CPM along with the diffusion and advection scheme. Section 3 
discusses various cases, including diffusion from two point sources, a continuous source and 
a moving source with either reflecting or absorbing boundary conditions and the flow profile 
in Poisueille's flow. Section 4, outlines future directions in developing the ADE scheme. We 
will address in a future paper certain additional conditions, flow with inertia, effects of low 
or high Peclet numbers, constitutive properties of the fluid phase, etc. As we have stated 
before, ours method provides flexibility and efficiency in the biologically relevant regime of 
low Peclet and Reynolds numbers and where high numerical accuracy is not crucial. 
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II. THE MODEL 



The Potts model is an energy-based, lattice Cellular-Automaton (CA) model equivalent 
to an Ising model with more than two degenerate spin values. We typically use a cubic 
lattice with periodic or fixed boundary condition in each direction. We use 3rd or 4th near- 
est neighbor interactions to reduce lattice-anisotropy induced alignment and pinning. Each 
lattice site in the Potts model has a spin value. The energy, or Hamiltonian sums the inter- 
action energies of these spins, according to predefined rules. In single-spin dynamics (like 
Metropolis dynamics) the spin lattice evolves towards its equilibrium state by minimizing 
the interaction energy through spin flips. Multi-spin dynamics like Kawasaki dynamics are 
also possible. Though the original Potts model studies focuses on equilibrium properties, 
it can also model quasi-equilibrium dynamical properties [38(. Using deterministic schemes 
for spin flips, patterns often stick in local minima. Finite-temperature Monte-Carlo schemes 
circumvent this problem. These schemes accept a spin flip with temperature dependent 
Boltzmann or modified Boltzmann probability if the co nfig uration encounters a potential 
barrier (a greater energy after the spin flip than before) 38] . 

We pick a target lattice site at random and one of its alien neighbor, also selected at 
random and attempt to flip the target spin to the value of the selected neighbor. In the 
modified Metropolis algorithm we employ, if the spin flip would produce a change in energy 
AH, we accept the change with probability P given by 

s f exp(-AH/T) if AH > 0, 
P{AH) = I y 1 ' (2) 
I 1 otherwise, 

where T is the fluctuation temperature. T controls the rate of acceptance of the proposed 
move. For very large T, all the moves are accepted and the dynamics is a random walk in the 
absence of barriers, i. e. interaction energies included in the Hamiltonian are effectively zero 
producing a disordered phase. For very small values of T the dynamics is deterministic and 
can trap in local minima. We choose T as the median value of the distribution of AH, which 
is below the order- disorder phase transition temperature. All of our results are very robust 
with respect to variation of T. S. Wong has recently shown that optimizing the dynamics of 
;he modified Metropolis algorithm requires changes to the acceptance probability in eqn. El 
39} ] . However, we do not implement these changes here. Our unit of time is Monte-Carlo 
sweep (MCS), where 1MCS = L 3 spin flip attempts, L being the system size. 



The CPM adapts the Potts model to the context of biology. A CPM cell is a collection 
of lattice sites with same spin value (or index) cxj. Each cell has a unique spin a (see fig. [I]). 
Cells may also have additional characteristics, e.g., a type r. Links between different sites 
with spins define cell boundaries. So cells have both volume and a surface area. The volume, 
area, radius relation is highly non-Euclidean for small cells. 
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FIG. 1: A typical cell configuration in CPM. The bold lines denote cell boundary. 

The CPM Hamiltonian contains a variable number of terms. The interaction between 
pairs of biological cells involves an adhesive or repulsive interfacial energy. This interfacial 
energy is precisely the Potts energy, the sum of the interactions of neighboring unlike spins, 
across a link. Each mismatched link contributes a cell-type dependent binding energy per 
unit area J(r, r ), where r and r are the type of cells on either side of the link. In the CPM 
the effective cell-cell interaction energy is: 

-^adhesion ^ J{j~cr(i ,j,k) j T<r(l,m,n) ) ( 1 ^(a(i,j,k),a(l,m,n)) ) j (3) 

(i ,j,k),(l,m,n)neighbors 

where 8, a i x = 1 if a = er, otherwise 8,^ s = 0. 

At any time t, a cell, of type, r, has a volume v(a,t) and surface area s(a,t). The vol- 
ume is simple to define, i>(<7 ) = ^ j k^(<ro,<r(i,j,k))i whereas surface area is more complex, 
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since it depends on the interaction range of the lattice, s(a ) = k ^(<y ,<j i3 ■ k ) f fc'(l ~ 
3(a(i,j,k),a(i' ,j' ,k')))i where (i,j,k) and (i',j',k') are neighboring lattice sites. Each cell has an 
effective volume elasticity, A„ and target volume v target (a,t). Larger values of A„, produce 
less compressible cells. We typically choose v 2 arget X v > other constraint energies. This com- 
pressibility makes little difference in low Reynolds number flow but makes pattern evolution 
less stiff. We also define a membrane elasticity, A s , and a target surface area s ta rget{.cr, t) to 
maintain the generalized shape of the cells. The energy contributions due to surface and 
volume fluctuations are: 

E S urface = \ s (s((T,t) — S targe t((J, t)) 2 , (4) 



Evolume = \ v {v(cr,t) - V target (<7, t)) 2 . (5) 

We can extend the Hamiltonian to include a uniform external force F, acting on all cells 
by including the term: 

Eforce ^ ^ F ' ^a(i,j,k),a(l,m,n)) i (6) 

(ijk),(lmn)neighbors 

where is the position vector at the lattice site (i, j, k). 

Previous CPM applications often treated fluid or ECM as a single large cell with no 
constraints. Here we take coarse-grained approach to describe ECM. We assume the ECM 
consists of hypothetical fluid cells (which we call particles to avoid confusion with the mod- 
eling of biological cells) having all the characteristic interactions and constraints of regular 
CPM cells. The volume constraint and the surface tension determine the elastic nature of the 
fluid. The fluid particles can move with respect to one another like regular CPM cell via spin 
flips. Thus, local pressure developed due to movement or enlargement of actual biological 
cells will translate into motion of the surrounding ECM. This fluid motion causes advection 
and mixing along with molecular diffusion of chemical morphogens. We restrict consider- 
ation to the highly overdamped viscous world that most biological cells experience, so our 
fluid particles lack inertia. We also restrict to situations where the velocity of movement is 
much less than the velocity of sound, which is one lattice unit per MCS. 

We now introduce a relative velocity constraint between the cells/particles which faith- 
fully captures the effects of shear due to the viscosity of the medium. In the CPM, velocity 
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is a cell property defined as the displacement of the center of mass of the cell per MCS. Since 
the velocity gradient terms in the direction of the flow (e.g., |^-) are the rate of change of 
volume [41| which eqn. 03 already includes, we need to keep only the contributions of cross 
terms of the form (ff 1 ) 2 - In an incompressible fluid (V ■ u = 0) the cross terms are the 
dissipation energy per unit volume p. In the CPM we model this term as : 



f -\ v-y-c (Vi*-Vj*) 2 \ (y l -y ] ) 2 + (z l -z J f 

^viscosity ^viscosity / J / J &%j ^2 \ 

i j l 3 V l 3 

+ cyclic permutation of (x, y, z), (7) 
where the j's are the indices of cells neighboring the zth cell and dij = 



a/ (xi — Xj) 2 + (yi — yj) 2 + (zi — Zj) 2 is the distance between the centers of cell i and cell 
j. Vi x is the x-component of the velocity of the ith cell. Since the cells are of irregular 
shape, we further weight the energy penalty by the cells common contact area Sij. We 
ensure that the cells are simply connected by using local connectivity checks during spin flip 
attempts. X V i SCO sity corresponds to the viscosity coefficient 77 in the Navier-Stokes equations. 
rj has dependence on other system parameters like J, X vo iume an d ^surface- 
The net Hamiltonian including fluids is then : 

H ^adhesion E sur j ace -\- Evolume ^viscosity (8) 



Diffusion Scheme 



Since the motion of the fluid particles takes care of advection, we need only to solve the 
diffusion on the current fluid particles configuration; _ D'S/ 2 C(x,t), where we have 

assumed that the diffusion constant D is constant and isotropic. Including an anisotropic 
or spatially varying D is a trivial extension of our method. We assume that the fluid 
particles carry chemical morphogens, whose distribution is uniform over a given fluid particle. 
Equivalently we can associate the chemical concentration with the center of mass of the 
particles and think of the diffusion as taking place between them (a FE view). Because the 
shape and volume of the fluid particles are irregular, their centers of mass do not correspond 
to lattice points. We then need a numerical scheme for diffusion among fluid particles. Few 
existing algorithms solve diffusion on a random or irregular lattice. The more sophisticated 
and accurate ones are computationally expensive j^j]. We use a naive iterated Euler method, 
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which is fast and stable and reproduces different biological experiments with good qualitative 
and fair quantitative accuracy. We could of course, use a more elaborate scheme, if necessary. 
We locally average the concentration among the particles nearest neighbors. The particle 
neighbors change as the fluid flows. We approximate the Laplacian V 2 C(r, t) (where C(r, t) 
is the chemical concentration and r is the center of mass coordinate of a fluid particle) by 



where R is the average radius of the fluid particles assuming them to be spherical (R = 
=1 N i — ) and Cj(r, t) is the concentration in nearest neighbor fluid particles. We use R 2 
instead of \ fl — rj| 2 to avoid Neumann instability. We can update the concentration once or 
multiple times per MCS depending on the diffusion coefficient of the chemical morphogen. 
D along with number of concentration updates per MCS controls the diffusion constant. For 
larger diffusion coefficients we update the diffusion step multiple times per MCS. For typical 
chemical morphogens like cAMP the diffusion constants are ~ 10~ 4 /im 2 /s , corresponding 
to a spread of 3-\/2.xl0~ 4 /z m (0.04/xm) per MCS. The typical lattice spacing in our CPM 
corresponds to 0.1 — 1.0/xm, hence a fluid particle has a width of ~ 0.3 — 3.0//m. Therefore 
a very small amount chemical needs to diffuse among adjacent fluid particle per MCS. In 
fact, for most practical purposes we need at most one diffusion step per five to ten MCS. 
Both reflecting and constant concentration (zero concentration is absorbing) boundaries are 
simple to implement. For reflecting boundary condition we exclude from the exchange of 
concentration (eqn. EJ any particle whose boundary is reflecting. For absorbing boundary 
condition (constant concentration), we define the particles whose boundaries are absorbing 
to always have zero concentration. The algorithm handles moving boundaries automatically 
and conserves chemical concentration, unlike many other methods. 

III. RESULTS AND DISCUSSION 
A. Poiseuille Flow 

We first discuss simulating Poiseuille flow of a viscous fluid in a cylindrical tube with rigid 
walls under a uniform force field (body force) using the CPM. We consider a cylindrical 
tube with a circular cross-section of radius 18 lattice units and length 200 lattice units 






i next to j 
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with periodic boundary condition along the length. Each fluid particle has an average 
volume of 64 lattice points. Our CPM AD scheme works for small forces when the fluid 
particles remain simply connected. Fig. |2] plots the ensemble averaged (50 different initial 
configurations) cross-sectional profile of the viscous flow at x — 150 and t = 2000 MCS 
with a small force F = 0.05x applied on all the fluid particles. The flow is in the x-direction 
only. The profile is parabolic as expected. We also show its fit to the analytical solution 




Y at Z = 20 



FIG. 2: (A) Velocity profile in a cylindrical flow under gravity and its fit to the analytical solution. 
(B) Velocity profile across a diameter of the tube with error bars, the bold segment along the x 
direction denotes the width of a fluid particle. 




FIG. 3: (A) The absolute error between the analytical solution V a naiyticai an d the CPM simulation 
^simulation at x = 150. (B) The rms error of over 50 ensembles. 

V(r) = ^R 2 {1 — (-^) 2 ), where rj is viscosity, is a fitting parameter, R — 18 the radius 
of the cylinder and r is the distance from the axis of the cylinder. What is surprising is 
the excellent agreement (an error of < 5% with the analytical result despite the coarseness 
of the simulation (only 9 particles across the diameter). All lattice points inside a fluid 
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particle have the same velocity, so regions where the shapes of the fluid particles are regular 
for most ensembles will produce plateaux in the velocity profile. Since in the CPM, the 
velocity of fluid particles is the center of mass velocity, a lattice point touching a boundary 
wall will have a small nonzero velocity as the center of mass of the corresponding particle 
lies in the interior (slip boundary). Since the energy contribution of the constraints near 
a boundary wall dominates the external applied force and the Monte-Carlo temperature, 
boundary particles are more regular in shape than interior particles. This regularity holds 
in all ensembles, hence the rms error in the velocity near a boundary wall is small as shown 
in fig. and the velocity just near the boundary in fig. El has a small plateau. Reducing 
the size of the fluid particless compared to the typical length scales of the flow reduces these 
anomalies. Increasing X V i SC osity increases the viscosity coefficient r\. A future paper will study 
the relation between rj and A, Monte-Carlo temperature, fluid particle size, etc. We shall 
also show additional biologically relevant hydrodynamic flows. 

We next verify our diffusion scheme under various biologically relevant conditions. CPM 
models using the modified finite temperature Metropolis algorithm have diffusion due to the 
movement of the CPM cells themselves which adds to the diffusion of chemical morphogens. 
However for most CPM simulations, e.g. a temperature 0.1 and other parameter values used 
through out this paper, the diffusion coefficient of the CPM cells is ~ 10~ A pixel 2 / MC S , 
which is much smaller than the typical diffusion coefficient of ~ 0.1, that we treat in this 
paper. Hence for pure diffusion in a static medium the fluid is effectively fixed. Besides 
the concentration profile of the diffusing chemical in a medium, diffusion in the presence of 
boundaries and moving source is crucially important in biology. We show that the results 
of our simple CPM diffusion from CPM agree very well with corresponding analytical cal- 
culations or finite element results. We also briefly discuss the validity of our method for 
diffusion in Poiseuille flow. 

The four cases we discuss below employ fluid particles with a target volume v tar get = 27, 
unless we mention otherwise. The lattice has 100x100x100 sites. We equilibrate and quench 
to remove any disconnected cells which the finite-temperature equilibration produces j^. 
In our model a chemical source at a given lattice point gives all lattice points which belong 
to the fluid particle containing that chosen point the same concentration. We apply one 
diffusion step per MCS. The first three cases are for static fluids. 
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B. Two Sources with Reflecting and Absorbing Boundaries 



We consider two point sources at 15, 50, 50 and 50, 50, 50 with initial concentra- 
tions (at t — 0) of 5 and 10 respectively. The chemicals diffuse from these in- 
stantaneous sources. The bounding planes of the cube are reflecting. Figure |U 
plots the corresponding one- dimensional diffusion profile projection (with no ensem- 
ble average) after elapsed times t = 50 MCS, 100 MCS, 200 MCS and fitted 

-(i-il) 2 -(2h\-x-x\) 2 -(2L2-x-xl) 2 

to the exact solution C(x,t) = ( ex P iDt + exp «Jt +exp &ot ) + 

,„ -(x-x2) 2 -(2L\-x-x2) 2 -(2L2-x-x2) 2 

y/~(4 Dt) ( ex P ADt "^^P ADt +exp 4Dt )■, D being a fit parameter. Here 
LI and L2 are the coordinates of the two reflecting boundaries and xl and x2 are the co- 
ordinates of the instantaneous sources. The diffusion profile matches matches very well at 
all times, even near the boundaries. The maximum relative error (compared to the exact 
solution) at t = 50MCS is 5%, again, surprisingly good for such a coarse simulation. We 
discuss errors in detail for absorbing boundary conditions in the next paragraph. The inset 
shows the variation of diffusion coefficient calculated from fitting to the exact solution as 
function of time. Though the diffusion constant should remain constant with time, we ob- 
serve a 5% — 6% variation from the asymptotic value at small times due to the sharp initial 
distribution producing structures smaller than the fluid particle scale of our coarse-grained 
scheme. We also studied the variation of the diffusion coefficient D as function of the target 
volume of the fluid particles. D is constant for variations over one order of magnitude of 
the target volume of the fluid particles as shown in the inset to fig. 03 Fig. El also shows 
the chemical profile for an absorbing barrier at x = and x = 100. The figures show that 
the diffusion coefficient for both reflecting and absorbing barriers cases is same. Hence for a 
wide range of fluid particle volumes, our CPM ADE algorithm faithfully reproduces diffusion 
from two point sources. 

Fig. IHlplots the relative error as a function of position at different times. As mentioned in 
the last paragraph the sharp distribution at the initial time produce errors large compared 
to later time, i.e., if we ignore the large errors in the tail (as the magnitude of concentration 
is extremely small in the tails) the relative error at t = 150 is < 5% whereas at t = 500 
it is < 0.5%. In the actual biological situation, cells secrete chemical morphogens over 
their whole membrane surface and hence such singular cases of high point concentration of 
chemical rarely occur. 
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FIG. 4: Symbols denote chemical profile projected on the rc-axis (summing in the Y and Z direc- 
tions) at t = 50 MCS, 100 MCS, 200 MCS for reflecting barriers at x = and 100 denoted by 
symbols. The solid line is a fit to the exact solution using the fitting parameter D. The inset shows 
variations in D with time, due to coarse graining and the initial sharp distribution. 

C. Two Sources with a Reflecting Obstacle inside the Medium 

Since biological cells can be impermeable to many chemical morphogens, they can act 
as reflecting boundaries within the fluid medium. We check this situation for the simple 
test case of a cylindrical barrier with rectangular cross-section (15x15, pixels centered at 
(50, 50, 50) and the axis along the x-direction) within the fluid medium. As in the previous 
case, we place two point sources near the two opposite corners of the rectangle at (50, 40, 40) 
and (50,60,60). We recover the same diffusion coefficient as in the previous case. Fig. [7| 
plots the one- dimensional cross-section of the diffusion profile at three different z positions. 
Along z = 40 and 60 the barrier is absent and the diffusion is merely the superposition of 
that from the two sources. Along z = 50 we see the effect of the reflecting barrier. The inset 
on the right side shows the concentration profile obtained from a finite element calculation 
for this situation which matches very well with the CPM concentration profile in the left 
inset. 
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FIG. 5: Chemical profile projected on the x-axis from diffusion of two instantaneous point sources 
for an absorbing barrier at x = and x = 100 for t = 70 MCS , 100 MCS , 200 MCS . The inset 
shows the variation of D with fluid particle volume on a log scale. 
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FIG. 6: Relative error along the x lattice direction for an absorbing barrier. 
D. Moving Continuous Source 

Moving cells often secrete morphogens. Hence we correctly simulate cells' chemotactic 
response to secreted chemicals only if we faithfully reproduce diffusion from moving sources. 
To test our simulation we assign an arbitrary fluid particle a constant concentration Co 
(continuous source) and uniform velocity v along the x direction. We keep the source 
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FIG. 7: Chemical profile cross-section along z = 40,50 and 60 in the presence of a cylindrical 
reflecting barrier with axis along the x direction and rectangular cross-section. The sources with 
initial concentration 5 and 10 are placed at 50, 40, 40 and 50, 60, 60. The inset shows 2d projection 
of the profile. The lower figure shows a two-dimensional projection obtained from finite element 
calculations. The concentration decreases as we go away from the center contours (yellow lines). 

sufficiently distant from the boundary to avoid boundary effects. At t = the source is 
at x = 35 and moves with velocity u = 0.04 pixel/MCS. We fit the CPM chemical profile 
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projection in the x direction (no ensemble average) with the ID analytic solution: 



a(x) = ((x-x ) 2 )/(AD), (10) 
(3 = u 2 /(AD), (11) 
j(x) = exp{(x-x )u/(2D)), (12) 




FIG. 8: CPM simulation of the chemical profile from a moving source (the other two coordinates 
integrated out) at t = 300 MCS. The solid line denotes the fit to eqn. ^1 . The inset shows a two 
dimensional projection of the simulation, where red denotes the highest chemical concentration 
and black the lowest. 

Figure |H1 shows that the CPM diffusion agrees very well with the analytical solution. The 
diffusion coefficient obtained from the fit is D = 0.29. 

E. Taylor-Aris Dispersion in Poiseuille flow 

We check the qualitative agreement of the dispersion coefficient obtained using our CPM 
ADE solver for Poiseuille flow along x direction (described in section 1111 AJ) in a cylindrical 
geometry. We compare our result for an initial delta distribution of chemical in the middle 
of the tube. After an initial transient, so that the chemical reaches the boundary in the 
y and z direction, we compare the diffusion coefficient of chemical distribution along the 
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x direction with the analytical result for the effective diffusion coefficient in Poiseuille flow 
along the cylinder axis, i.e. D effective = -Do + f^§~- Here D is the diffusion coefficient in the 
absence of flow, v is the average velocity and R is the radius of the cylinder. Figure El shows 
our results where we have plotted the variation of effective diffusion coefficient with time (in 
MCS) for different mean velocities. For the analytical results given in solid lines in fig. 
we use the mean velocity obtained from the fit of the velocity profile like as shown in fig. |21 
Figure ITU1 shows snapshots of chemical profile after 500 MCS for D Q = and D = 0.12. 
In this case, we start with a thin sheet of chemical in the x — z plane, at x = 75, i.e., 
all the particles at x = 75 have uniform concentration. We take a cut at y — 20 to see 
the evolved chemical profile. A detailed study of evolution of chemical profile, effect of 
embedded objects like sphere under both stationary and moving condition will be reported 
in our future communication. 
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FIG. 9: Effective diffusion coefficient for mean flow velocities v = 0.056, 0.028 and 0.016 respec- 
tively from top to the bottom curve. The solid lines show the analytical results corresponding to 
above mean velocities. 



IV. CONCLUSIONS 

We have implemented fluid flow, advection and diffusion in the framework of the CPM, 
avoiding the programming complexity and computational demands associated with imple- 
menting a finite-element or finite-difference Navier-Stokes simulation and interfacing it with 
the CPM lattice. 
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FIG. 10: Cross-sectional profile (x — z plane, y = 20) of chemical concentration at t = 500MCS 
starting from an initial set of particles at x = 75 with uniform concentration; (A) with no diffusion, 
(B) with diffusion. Particles which move away from the plane in the y direction cannot be seen in 
the figures. 

We have used three biologically relevant test cases to verify our method. All our results for 
diffusion in the presence of boundaries or moving sources agree very well with corresponding 
analytical or finite-element solutions. The errors in our scheme are large if we try to probe 
far below the diffusion time scale or the fluid particle length scale, but the results are 
qualitatively correct. Thus we must be cautious when applying this scheme to large Pe 
number flows. The requirement that fluid particles remains connected, limits the method 
to low Re. However, since most biological mechanisms operate at low Re our CPM ADE 
solver is appropriate for many cell-level simulations. 
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